Probing Majorana edge states with a flux qubit 
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A pair of counter-propagating Majorana edge modes appears in chiral p-wave superconductors 
and in other superconducting systems belonging to the same universahty class. These modes can be 
described by an Ising conformal field theory. We show how a superconducting flux qubit attached 
to such a system couples to the two chiral edge modes via the disorder field of the Ising model. 
Due to this coupling, measuring the back-action of the edge states on the qubit allows to probe the 
properties of Majorana edge modes. 
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I. INTRODUCTION 

Chiral Majorana fermion edge states were originally 
predicted to exist in the 5/2 fractional quantum Hall 
plateauii These edge states support not only neutral 
fermionic excitations but also more exotic edge vortices. 
A single edge vortex corresponds to a tt phase shift to 
all fermions situated to one side of it<^"— Two edge vor- 
tices may cither fuse into an edge fermion or annihilate 
each other, with the outcome depending on the preceding 
evolution of the system. In other words, the edge theory 
(together with the corresponding bulk theory) possesses 
non-Abelian statistics!^"— This unusual physics and its 
potential applications to topological quantum computa- 
tion are the reasons why the Majorana edge states have 
attracted much attention recentlyi^"— 

Similar non-Abelian anyons and their corresponding 
edge states appear in superconducting systems as well. 
Initially it was discovered that p-wave superconductors 
support non-Abelian anyons in the bulk and chiral Ma- 
jorana edge states.— li^i^ Later it was shown that de- 
positing a conventional s-wave superconductor on the 
surface of a topological insulator while breaking time- 
reversal symmetry provides an alternative route to real- 
ize these non-Abelian states.—"— Alternative proposals 
include substituting the topological insulator by a two- 
dimensional electron gas with spin-orbit coupling22r— or 
by a half-metal.— The realizations of Majorana edge 
states using s-wave superconductors have the following 
advantages: first, they rely on combining simple, well- 
studied ingredients. Second, the materials do not have to 
be extremely pure unlike samples needed to support the 
fractional quantum Hall edge states. Finally, the super- 
conducting implementations of Majorana fermions may 
feature a larger bulk excitation gap and may therefore be 
operated at higher temperatures. 

The downside of the superconducting implementations 
of Majorana edge states is the lack of means to manipu- 
late edge vortices >i2ii^ Different from the 5/2 fractional 
quantum Hall state, the edge vortices are not coupled 
to charge and thus cannot be controlled by applying 
voltagesi^ Therefore, the standard proposal to probe 



the edge vortices in superconducting systems is to in- 
ject fermion excitations into the edge, to let them split 
into edge vortices, and finally to conclude about the be- 
havior of the edge vortices from the detection of the 
fermion excitations after the subsequent fusion of edge 
vorticesM^i^e 

In this paper, we propose a more direct way to manipu- 
late and measure edge vortices using a flux qubit consist- 
ing of a superconducting ring interrupted by a Josephson 
junction.— Our main idea is based on the following ob- 
servations: first, an edge vortex is created when a super- 
conducting vortex crosses the edge. Second, the motion 
of the superconducting vortices can be fully controlled by 
a flux qubit, since by applying a flux bias to the qubit 
one can tune the energy cost for a vortex being present 
in the superconducting ringj^ In this way, attaching a 
fiux qubit to a system supporting Majorana edge states 
allows to directly create, control, and measure edge vor- 
tices without relying on splitting and fusing fermionic 
excitations. 

We note that our proposal is not necessarily advanta- 
geous for the purposes of topological quantum computing 
since quantum computing with Majorana fermions may 
even be realized without ever using edge statesi^"— In- 
stead the aim of our investigation is to develop a better 
tool for probing the fractional excitations of the edge the- 
ory. 

The paper is organized as follows. In Sec. |TT1 we dis- 
cuss a schematic setup of a system where a pair of chiral 
Majorana fermion edge modes couple to a fiux qubit as 
a probe of the edge states and briefly list our main find- 
ings. In Sec. mil we review the connection between the 
one-dimensional critical transverse-fleld Ising model and 
Majorana fermion modes. We identify the vortex tun- 
neling operators between two edge states as the disorder 
fields of the Ising model, and subsequently derive an effec- 
tive Hamiltonian for the fiux qubit coupled to Majorana 
modes. In Sec. IIV| we provide the necessary formalism 
for evaluating the expectation values for the flux qubit 
state and qubit susceptibilities. In Sec. |V] and Sec. I VII 
we compute the qubit expectation values and the two- 
point qubit correlation functions in the presence of the 
edge state coupling, and use these results to derive the 



qubit susceptibility. In Sec. I VIII we analyze higher or- 
der corrections to correlation functions of the qubit state. 
We summarize our results in Sec. IVIIII Additionally, we 
provide a brief overview of the flux qubit Hamiltonian 
in App. |X1 In App. |B] we reduce the flux qubit Hamil- 
tonian to that of a two-level system and derive the cou- 
pling between the flux qubit and the Majorana modes. 
In App. [Cl we give the form of the four point correla- 
tion function for the disorder field of the Ising model. 
Finally in App. [Dl we provide the detailed derivation of 
the higher order corrections to the correlation functions 
of the qubit state. 



flux qubit 




II. SETUP OF THE SYSTEM 

In this work, wc consider the following setup: a strip 
of s-wave superconductor is deposited on the surface of 
either a three-dimensional topological insulator or a semi- 
conductor with strong spin-orbit coupling and broken 
time-reversal symmetry (or any other superconducting 
setup supporting Majorana edge states). As depicted in 
Fig. [U a pair of counter-propagating Majorana fermion 
edge modes appears at the two opposite edges of the 
superconductor jiSii^ To avoid mixing between counter- 
propagating edge states, the width of the superconduc- 
tor should be much larger than the superconducting co- 
herence length hvp/^- Here and in the following, vp 
denotes the Fermi velocity of the topological insulator 
(semiconductor) and A the proximity-induced supercon- 
ducting pair-potential. In order to avoid mixing of the 
two counter-propagating edge modes at the ends of the 
sample, we require the length of the superconducting 
strip to be longer than the dephasing length. 

A flux qubit, consisting of a superconducting ring with 
a small inductance interrupted by a Josephson junction, 
is attached to the heterostructure supporting the Majo- 
rana edge modes, as shown in Fig. [TJ By applying an 
external flux $, the two classical states of the supercon- 
ducting ring corresponding to the phase difference of 
and 2tt across the junction can be tuned to be almost de- 
generate.— In this regime, the flux qubit can be viewed 
as a quantum two-level system with an energy difference 
e (which we choose to be positive) between the states |0) 
and |27r) and a tunneling amplitude 5 between them. As 
described in App. [XI the energy difference e can be easily 
tuned by the external flux $ threaded through the ring. 

The transition between the two qubit states is equiv- 
alent to the process of a vortex tunneling through the 
Josephson junction in or out of the superconducting 
ring. For convenience, we will refer to the Hilbert space 
spanned by the qubit states |0) and |27r) as a spin-I/2 
system. For example, we are going to call the Pauli ma- 
trices CT^'^'^ acting on the qubit states the qubit spin. 

A vortex tunneling through the weak link in the super- 
conductor from one edge to the other is a phase slip of 2-k 
of the superconducting phase difference at the tunneling 
point. Due to this event, all fermions to one side of the 



FIG. 1: Schematic setup of the Majorana fermion edge modes 
coupled to a flux qubit. A pair of counter-propagating edge 
modes appears at two opposite edges of a topological super- 
conductor. A flux qubit, that consists of a superconducting 
ring and a Josephson junction, shown as a gray rectangle, is 
attached to the the superconductor in such a way that it does 
not interrupt the edge states flow. As indicated by the arrow 
across the weak link, vortices can tunnel in and out of the 
superconducting ring through the Josephson junction. 

weak link gain a phase of tt. As will be shown below, the 
vortex tunneling operator can be identified with the op- 
erator of the disorder field of a one-dimensional critical 
Ising model onto which the Majorana edge modes can be 
mapped. 

Since vortex tunneling events couple the qubit spin to 
the Majorana edge modes, we expect various observables 
of the qubit to carry signatures of this coupling. The 
main theory parameter that we are after is the scaling 
dimension = 1/8 of the edge vortex operator (disor- 
der field). Our main results apply to the regime when 
vortex tunneling is weak e ^ (5. 

We find that the reduction of the spin expectation 
value in the z-direction due to the vortex tunneling ac- 
quires a nontrivial scaling exponent 



Similarly, the spin expectation value along the x- 
direction is proportional to £2A,,-i _ ^-3/4 ^j-^ereby prob- 
ing the scaling dimension of the disorder field. 

The finite frequency susceptibilities that characterize 
the response of the polarization of the qubit spin to a 
perturbation with frequency lo provide additional infor- 
mation about the Majorana edge states. The suscepti- 
bility Xzz(w), which characterizes the change of (ctz) due 
to a modulation of with frequency lu, is measurable 
with current experimental techniques. It can be mea- 
sured by modulating the external flux $ and reading out 
the current from a dc-SQUID coupled to the qubit i^iH 

The frequency dependence of the susceptibilities ex- 
hibits a non-Lorentzian resonant response around the fre- 
quency w « e (here and in the following, we set h = \). It 
is modified by the coupling to the Majorana edge states 
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and shows the scaling behavior 
lx(w)| oc 



1-2A, 



|w-e|3/4' 



(2) 



as long as e ^ |w — ej, and the distance |a; — e\ from the 
resonance is larger than the width of the resonance. The 
phase change of susceptibility at the resonance Scf) = 3tt/A 
is different from the tt phase change for a usual oscillator. 
The origin of the extra tt/4 phase shift is the Abelian part 
of the statistical angle of the vortex excitations<^ 



III. EDGE STATES AND COUPLING TO THE 
QUBIT 

A. Coupling of the flux qubit to the edge states 

The flux qubit has two low energy states, correspond- 
ing to a phase difference = or </> = 27r across the 
Josephson junction at x = xq. The Hamiltonian of the 
qubit is given by 



Q 



(3) 



The energy difference e can be tuned by applying an ex- 
ternal flux to the qubit while the tunneling amplitude 
^ > can be manipulated by changing the Josephson 
coupling of the junction<21 As discussed in App. the 
two levels described in ([3]) represent the two lowest energy 
states localized at the two energy minima of a double well 
potential. In order for the two-level approximation to be 
accurate, the energies d, e as well as the driving frequency 
Lu have to be much smaller than the level spacing at each 
well. The tunneling phase a is proportional to the charge 
induced on the sides of the junction and its fluctuations 
are the main source of qubit decoherence. For simplicity 
we neglect the charge noise so that we can assume that 
a is static and set it to zero without loss of generality. 
The qubit Hamiltonian now reads 



Tin ~ <J CF 



(4) 



When there is no phase difference across the Josephson 
junction (</) = 0), the Hamiltonian of the chiral Majorana 
modes appearing at the edges of the superconductor, as 
shown in Fig. [1] reads 



where vm is the velocity of the Majorana modes, and 
ipuix) and ipd{x) are the Majorana fermion fields at the 
upper and lower edges of the superconductor in Fig. [T] 
The sign difference between the terms containing ip^ 
and is due to the fact that the modes are counter- 
propagating. The Majorana fermion fields obey the anti- 
commutation relations 



{^l)u{x),ipu{x')) ^ {i/jd{x),'ipd{x')} = 2TrS{x - x' 
{Mx),Mx')} = Q- 



(6) 



A vortex tunneling through the weak link at x = Xq 
advances the phase of each Cooper pair in the region 
X < Xq hy 27r. For Majorana fermions, just like any 
other fermions, this results in phase shift of tt. The effect 
of this phase shift is a gauge transformation 



where the parity operator P is given by 



P = exp 



ITT / dx Pe{x) 



(7) 



(8) 



with the fermion density Pe{x) = ijj^ {x)ip{x) and ip = 
(jpu + iipd) /2y/7r. We refer to App. |B] for a derivation 
of the qubit Hamiltonian and the gauge transformation 
Eq. ([7]). When the phase difference between two sides of 
the Josephson junction is exactly tt, the Majorana modes 
approaching the junction are fully reflectedJ^ Since this 
phenomenon occurs only very close to the phase differ- 
ence of TT, where the system only spends a short amount 
of time during the process of a phase slip, we will ne- 
glect the effect of this backscattering. The relation be- 
tween the phase slip and the parity operator was dis- 
cussed and used in previous work focusing on the 5/2 
fractional quantum Hall statei^i^ii^ 

Combining the Hamiltonian of the Majorana edge 
states ([51 [7|) with the qubit Hamiltonian (|4]) , we get the 
full Hamiltonian of the coupled system in the basis of |0) 
and |27r): 



n 



Hmf 
PHmfP 



He 



(9) 



ipuix)dx'ipuix)], (5) pointi 



The first part of Hamiltonian represents the chiral Ma- 
jorana edge states coupled to the phase slip of the su- 
perconductor while the second part is the bare flux qubit 
Hamiltonian. 

Because the parity operator ([5]) is highly nonlocal if 
expressed in terms of Majorana fermions, it is desirable 
to map the Majorana modes on a system where the vor- 
tex tunneling event becomes a local operator. To this 
end, we establish the equivalence of the chiral Majorana 
edge modes with the long wavelength limit of the one- 
dimensional transverse-field Ising model at its critical 
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B. Mapping on the critical Ising model 

The lattice Hamiltonian of the Ising model at the crit- 
ical point is given by^i^ 



Hr 



(10) 
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where s" are the spin-1/2 operators at site n. With the 
Jordan- Wigner transformation, 



Vet 



=4 exp(i7r2]]c]cj) 



l-2cl, 



(11) 



j<n 



the Ising model pH)) can be cast in terms of fermions as 
Hi ^ J Y^lic^ - cl)(c„+i + c,\+i) + 4 c„ - c„4]. (12) 

n 

Here sf = (sf ± isf)/2 obey the usual onsite spin com- 
mutation relations while the fermions operators c\ and 
obey canonical anti-commutation relations. 

For each fermion, we introduce a pair of Majorana op- 
erators = tpj^ and i/'n = n such that 



,-i7r/4 



(13) 



The Majorana fermions satisfy the Clifford algebra 

{^m, Ipn} = {ipm, ^n} = 2J„i„, {'(/'m, Ipn} = 0. (14) 

In terms of the Majorana operators, the Hamiltonian 
assumes the form 



Hr 



iJ X ^ 

y ^ 



(V'n^/'n+l - ■0nV'«+l 

+ V'nV'n+1 - ■0nV'n-|-l " ^Ipni^n)- (15) 



In the long wavelength limit, the Hamiltonian (jlSp re- 
duces to (O with the identification of the continuum Ma- 
jorana operators 



■011 (x ) h-s- ./-?/!„, ipci{x) ^ \ -i'n, x^na (16) 
V a V <^ 

and the velocity vm 2Ja. To complete the map- 
ping, the bandwidth of the Ising model should be re- 
lated to the cutoff energy A of the linear dispersion of 
the Majorana edge states, A i— >■ J. Thereby, a pair of 
counter-propagating Majorana edge states, ipu{x) and 
ipd{x), can be mapped on the low energy sector of the 
one-dimensional transverse-field Ising model at its criti- 
cal point. 

For the parity operator ([5]) , we obtain a representation 
in terms of the Ising model with the following procedure: 
we first discretize / °dx pe{x) using the mapping 
and identify = no a as a lattice point on the Ising 
model. Thereafter, we obtain an expression for the vortex 
tunneling operator P in terms of the Ising model 

P ^ exp(i7r = n *l = <o+i/2' (17) 

i<"o i<"o 

by using Eq. and the Jordan- Wigner transformation 
PD) . Here, is the disorder field of the Ising model, i.e.. 



the dual field of the spin field<^"— The Ising Hamiltonian 
has a form identical to Eq. pU)) when expressed through 
H operators. 



-1/2 



-1/2) 



(18) 



with ~ ^n^n+i'^ We see that the parity operator 

is indeed a local operator in the dual description of the 
Ising model. After mapping on the Ising model Eq. ([7]) 
becomes (here and in the following, we use the shortcut 
notation fi = 



PHmfP ^ Mno + l/2 Hi /^,io + l/2, 



(19) 



and the full Hamiltonian of Majorana edge states and the 
fiux qubit ^ maps onto 



H^Hi^ 



Hr 







Ai,io + l/2 Hi Mno + l/2 

Finally, an additional unitary transformation 
Hi ^ VHiV\ 



Hr 



V ^V^ = 



1 

+ 



yields 



e , S ^ 
2 • 2''' ^"0+1/2- 



(20) 

(21) 
(22) 

(23) 



Here, are the Pauli matrices acting in the Hilbert space 
spanned by |0) and yUno+i/2|27i')- The operators of the 
qubit spin can be expressed through t^'V'^ as 



T>„o+i/2- (24) 



Mno + 1/2, 



We use the Hamiltonian in the form of Eq. ([23| and the 
qubit spin operators (|24p in the rest of the paper. 

The way of identifying two edge Majorana states with 
a complete transverse field Ising model presented above 
is different from the one commonly used in preceding re- 
search. Usually, the chiral part of the Ising model is iden- 
tified with a single Majorana edgci^^ The advantages of 
our method are the possibility to write a complete Hamil- 
tonian of the problem and simplified book-keeping, while 
its drawback is the need for the right-moving edge and 
the left-moving edge to have the same geometries. Over- 
all the differences are not important and both methods 
can be used interchangeably. 



IV. FORMALISM 

To probe the universal properties of Majorana edge 
states, the energy scales of the qubit should be much 
smaller than the cutoff scale of the Ising model, e, 6 A. 
In the weak coupling limit e ^ (5, we construct a per- 
turbation theory in S/e by separating the Hamiltonian 
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Hi = Hq + V into an unperturbed part and a perturba- 
tion 



Ho 



V 



(25) 



Without loss of generality we set e > 0, so that the 
ground state of the unperturbed qubit is |0). For brevity 
we omit the spatial coordinate of the fi operator in the 
following since it is always the same in the setup that we 
consider. 

We use the interaction picture with time-dependent 
operators 



o{t) = 



J'Hnt 



Oe 



-iHot 



The perturbation V{t) in this picture is given by 



{t)+T-{t% 



(26) 



(27) 



where t {t) = t are the time-dependent raising 
and lowering operators. The structure of the raising and 
lowering operators leads to physics similar to the Kondo 
and Luttinger liquid resonant tunneling problems j^i^i^ 
In the calculation we need the real-time two-point and 
four-point correlation functions of /i in the long-time limit 
A\t — t'\ ^ 1. The two-point correlation function is 



{^l{t)^,{t')) 



g-«sgn(t-t')7r/8 
A2A^|t-f'|2A, 



(28) 



where sgn(x) denotes the sign of x, and = 1/8 the 
scaling dimension of the n field.— The phase shift 7r/8 of 
the two-point correlator is the Abelian part of the sta- 
tistical angle for the Ising anyons braiding rules Corre- 
lation functions involving combination of multiple fields 
can be obtained via the underlying Ising conformal field 
theory or via a bosonization scheme42r— The expression 
for the four-point correlation function is given in App. [C] 
due to its length. For brevity we will measure energies 
in units of A and times in units of 1/A in the following 
calculation and restore the dimensionality in the final re- 
sult. 

We are interested in observables of the fiux qubit: 
the spin expectation values and the spin susceptibilities. 
We use time-dependent perturbation theory to calculate 
these quantities!^ This method is straightforward be- 
cause of the simple form of the perturbing Hamiltonian 
(|27p in terms of raising and lowering operators. 

Assuming that the system is in the unperturbed 
ground state at time to ^ — oo, the expectation value 
of a qubit spin operator cr"(t) is expressed through the 
S-matrix S(t, t'), 

{a"it)) = {S{t,to)^a"it)Sit,to))o, (29) 

S{t, t') = r exp (^-i V{s)ds^ ,t>t'. (30) 

Here, T is the time-ordering operator and (•)o is the ex- 
pectation value with respect to the unperturbed ground 



state. Similarly, the two-point correlation functions of 
the qubit spin are given by 

K(t)a'3(0)) = {S\t,to)<J^{t)S{t,0)a^{0)S{0,to))o. 

(31) 

The perturbative calculation for both the expectation 
values and correlation functions is done by expanding 
the 5'-matrices in V order by order. This procedure is 
equivalent to the Schwinger-Keldysh formalism with the 
expansion of S and S'^ corresponding to insertions on the 
forward and backward Keldysh contour. 

According to linear response theory, the susceptibility 
is given by the Fourier transform of the retarded corre- 
lation function of the qubit^ 

Xc.p{^)^q dte-*(K(t),a^(0)]), 

= -2/ dte'"*Im(CT"(t)CT^(0))c, (32) 
Jo 

where {■)c denotes the cumulant, 

K(t)a^(0))e = (a"(t)a^(0)) - (a"(t))(a^(0)), (33) 

and we have used (cr'5(0)cr"(i))c = (cr"(t)cr^(0))*. We 
see that in order to calculate the susceptibilities only the 
imaginary part of the correlation functions for t > is 
required. 



V. EXPECTATION VALUES OF THE QUBIT 
SPIN 

In this section, we calculate the expectation values of 
the qubit spin due to coupling with the Majorana edge 
states to the lowest non- vanishing order. Using the iden- 
tity 



1 - 2CT"cr+, 



we obtain 



(34) 



(35) 



since (ctz)'"' = 1. 

The first non- vanishing correction in the perturbative 
calculation of {a~a~^) is of second order in V. By ex- 
panding S and S'^ in Eq. (f^ . we obtain 



^0 

dh / d^2/^ (36) 

-oo J — oo 
r = {Vit2)T~T+V{tl))o. 

The integrand originates from the first order expansion 
of both 5* and S't . The second order contributions from 
the same S- or S't-matrix vanish due to the structure of 
V in the c^ubit spin space. 

Substituting (|27|) and (|28)) into the integrand yields 



r 



J2gi£(ii-t2)-isgn(t2-ti)7i 



(37) 
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By evaluating the integral in Eq. (|36p . we find 



8£ 



(38) 



where r(a::) denotes the Gamma function. 

The expectation value of in the unperturbed ground 
state vanishes. The first non-vanishing contribution to 
(cr^) arises to first order in 6/e. Expanding 5* and in 
Eq. ((^ to the first order yields 



r 



-z([r>(0),\/(ti)])o = 



dtiT, (39) 
sin(-rfi + f )J 



l^il 



2A„ 



after substituting cr^ from Eq. ([24|) and employing the 
two point correlator, Eq. ([25)1 . Evaluating we find 



(40) 



rl-2A„ 



Finally, (c^) = to all orders in perturbation theory 
since the Hamiltonian is invariant under i— > — ct^. 



VI. CORRELATION FUNCTIONS AND 
SUSCEPTIBILITIES OF THE FLUX QUBIT SPIN 

Since we are interested in the behavior of susceptibili- 
ties at frequencies close to the resonance w ~ e, we only 
need to obtain the long-time asymptotic of the correla- 
tion functions of the qubit spin. Using (|24|) and ((28)) . we 
immediately obtain that 



(a-(t)a-(O)), = 



2A„ 



(41) 



is non-vanishing to zeroth order. This is due to the 
fact that flipping the qubit spin automatically involves 
creation of an edge vortex, and is exactly the spin 
flip operator. In the same manner, one obtains that 
{ay{t)(Ty{Q))c = (CT^(t)cr^(0))c to zeroth order. 

Concentrating next on the mixed correlator, the rela- 
tions (IMl) and (IMl) yield 



(a-(i)a^(0))c = -2W)T%t)T-{Q)T+{Q)W (42) 

The leading non-vanishing term in this correlation func- 
tion is of flrst order in 5 and given by 



(43) 



in the long-time limit. 

The leading order contribution to ((T^(t)cr^(0))c can 
be evaluated using (j28p with expansions of S and to 
second order in 5. In the long-time limit, the leading 
contribution of the correlation function is given by 



(a^(0a^(0))(2) 



KW^"(0))c. 



(44) 



Correlators containing a single vanish because of the 
invariance under i— >■ —a^. We see that all the non- 
vanishing two-point correlation functions are the same 
up to overall prefactors. Therefore, we will focus on 
(cr^(t)CT^(0))c in the following. 



A. Energy renormalization and damping 

The coupling of the flux qubit to the continuum Majo- 
rana edge states can be thought of as a two-level system 
coupled to an environment via the interaction ([?7)) . This 
coupling leads to self-energy corrections E for the qubit 
Hamiltonian 



(45) 



that effectively shifts the energy spectrum and can also 
induce damping 4^ Since we are only interested in qubit 
observables, we focus on the structure of S for the two- 
level system and do not discuss the self-energy correction 
of the Majorana edge states. 

To second order, the self-energy correction for two spin 
states can be written in terms of the perturbed Hamilto- 



nian as^ 



= (a; 0|y + V{E^ + ^0+ - Ho)" V|0; /?), (46) 

where Ea is the energy for the spin-a =t, i qubit states 
and |a; 0) indicates that the Ising model is in its ground 
state with spin-a for the qubit state. Due to the structure 
of the Hamiltonian p5|) . the first order correction to the 
self-energy vanishes. Additionally, the off-diagonal self- 
energy corrections vanish also to second order. 

By inserting a complete set ^ \Ej] I3){f3; Ei\ — 1 
of the Hilbert space of Ho with Ej denoting the complete 
set of eigenstates with energy Ej for the Ising sector, the 
diagonal elements of the self-energy become 



y 



E 

Ei,l3 



{a;0\V\Ei;^){P;Ei\V\0;c 
Ea + iO+ - [Ei + Ep) 



(47) 



Because V — ~{d /2)TxfJ., only terms with a ^ (3 give 
non-vanishing contributions such that 



\f,\Ej){Ei\fx\0) 
Ei+iO+ 



(48) 



where -I- corresponds to a =1, and — to a =t- The 
diagonal elements of the self-energy in Eq. can be 
cast to the form 



(M(i)M(O)). 



(49) 



To see that (|49)) is equal to p8)) . we first insert a complete 
set of states of the Ising model, then write the time evo- 
lution of fj, in the Heisenberg picture, and finally evaluate 
the integral. 
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Evaluating Eq. gg]) with Eq. ^ yields 



■"tt 



(50) 



where we have used £ > 0. The absence of the imaginary 
part for E-j-t- indicates that the spin-up state is stable. 
The self-energy thus gives an energy shift to the spin-up 
state while it gives an energy shift with a damping to the 
spin-down state, 



e £ 
2 2 



(51) 



The energy renormalization and damping (j5ip alter the 
time evolution of the ground state correlation function 



(r+(t)r-(0))o = e' 



-i(e+v)t--it/2 



(52) 



where the energy renormalization and damping v — 
17/2 = Sj,! - are given by 



cos2(|)r(|)52 



2e 



1-2A„ 
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r(!)<52 



2V2£ 



1-2A„ 



(53) 



At zero temperature, this correlator is the only non- 
vanishing qubit correlator that enters in the perturbative 
calculation. Therefore, the effect of the self energy can 
be captured by replacing 



e ^ e - 



7, 



(54) 



in the qubit correlation functions computed in the long- 
time limit excluding the self-energy correction. Using the 
replacement rule ()54p . one obtains the zero temperature 
correlator 



(a-(t)a-(O)), 



-i{e+v)t — it/2-iTT/S 



t 



2A„ 



(55) 



The energy renormalization and the induced damping 
(j51|l do not arise explicitly in the lowest-order perturba- 
tion and require the resummation of the most divergent 
contributions to all orders in perturbation theory. In 
a system where Wick's theorem applies, the resumma- 
tion for the self-energy can be derived explicitly from a 
diagrammatic perturbation scheme4^ Because the corre- 
lation functions of multiple ^'s do not obey the Wick's 
theorem (see App. [C|), the resummation procedure for 
our system becomes more complicated. In the long time 
limit, however, the most divergent contributions in all 
orders can be collected by using the operator product 
expansion for two fields that resembles the structure of 
the Wick's tlieorem4ii^ 



B. Finite temperature 

Besides 7, finite temperature is an alternative source of 
decoherence. The finite temperature correlators of disor- 
der fields are readily obtained from the zero temperature 
correlators using a conformal transformatioii^ 



1 



1-^ 



[TTkBT] 



2A„ 



[sinh(7rfcBTt)]2A^ 



(56) 



where T denotes temperature and fc^ the Boltzmann con- 
stant. The finite temperature correlator (CT^(t)cr^(0))c 
in the long-time limit can be obtained by substituting 
Eq. ([55| into ([55]) with the proviso e ^ fcsT such that 
the temperature has no direct effect on the qubit dynam- 
ics. 



C. Susceptibility 

With the correlation functions derived above, we are 
now in the position to evaluate susceptibilities of the 
qubit. We should keep in mind that these correlators are 
valid only in the long-time limit and can only be used 
to study the behavior of the susceptibilities close to the 
resonant frequency cj « e. 

Evaluating Eq. ((32|) with Eq. ([55)1 yields the suscepti- 
bility at zero temperature around the resonance. 



„i3 



[i(e + t/-w)-h7/2]i- 



2A„ 



(57) 



where v and 7 are given in (|53p . Here, we note that the 
susceptibility (|57p shows non-Lorentzian response. This 
is in contrast to the conventional Lorentzian response of a 
two- level system weakly coupled to the environment 
If we neglect v and 7, which are of higher order in S/e, 
this susceptibility reduces to 



r(!) 



for uj < e, 



\UJ - £|1-2A„ I giS^/i^ for LU > 



(58) 



so it diverges and changes the phase by 37r/4 at the res- 
onant frequency. We can attribute this phase change to 
the phase shift of the correlator of two disorder fields in 
Eq. (USD. 

The presence of damping 7 in Eq. ([57)) provides a cutoff 
for the divergence of the response on resonance. The 
maximal susceptibility is reached at uj = e + v, and its 
value is given by 



\Xxxie + iy)\ = 



2l-2A,p(3) 
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1-2A„ 



(59) 



Using the proportionality of the correlation functions 
(gSl) and gl]), one gets that Xxz = Xzx = -{5/e)xxx 
and Xzz = iS/£)'^Xxx- It is interesting to note that when 
(5^0 both 

Xxx and Xxz are divergent while Xzz vanishes 
at the resonance. 

In Fig. [21 the absolute value of the susceptibility 
|Xa:a;(<^)l closc to thc rcsonance is plotted as a function of 
frequency. The dotted line shows the modulus of Eq. ([58[) 
for 1/ = 7 = while the dashed line shows that of 
Eq. ()57p . A renormalization of the resonant frequency 
v becomes clearly visible when comparing the peak posi- 
tion of the dashed line to that of the dotted line. 

The conformal dimension of the vortex excitation can 
be measured in thc region with e ^ |a; — ej > 7 where 



\Xxx{,Uj)\ 



r(|) 



|1-2A„ 



(60) 
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FIG. 2: Plot of the magnitude of the susceptibility jxa:a:('^)l 
as a function of frequency cu close to resonance e. The dotted 
line shows the zero temperature susceptibility in the absence 
of the damping and energy renormalization while the dashed 
line shows the result in the presence of the energy shift and 
the damping in Eq. (|5H) . The parameters used for the plot 
are e — O.IA and 5/e = 0.2. The solid line shows a plot of 
the finite temperature susceptibility with knT = 0.02e. 



As discussed in Appendix [Cl the four-point correlator, 
in principle, contains information about the non-Abelian 
statistics of the particles because changing the order of 
the fields in the correlation function not only alters the 
phase but can also change the functional form of the cor- 
relatorfi^ It is thus interesting to go beyond the low- 
est non-vanishing order. Additionally, doing so allows to 
check the consistency of the calculation of the self-energy 
correction done in Sec. IVI Al 

As an example we focus on the second order correc- 
tion to the ((T^ {t)<J^ (0))c correlator in the long-time limit. 
The details of the calculation are given in App. |D] and 
the result in Eq. (|D33p . The dominant correction is a 
power law divergence 



l-{ii. + l)t , (62) 



which is just the second order in 5 expansion of the mod- 
ified correlation function 



-i{e+u)t --ft/2 -iTr/8 

(a-(t)a-(O)), (X . (63) 



t 



Moreover, both Xxz and Xzz exhibit the same scaling 
behavior. 

The finite temperature susceptibility of Xxxi^^iT) can 
be evaluated from the correlation function ([55]) subjected 
to the transformation ([55]) . The result is plotted as the 
solid line in Fig. [5] An immediate effect of the tempera- 
ture is that it also introduces a cutoff for the divergence 
on resonance. For instance, the resonance peak of the 
susceptibility yields a different scaling behavior with re- 
spect to the temperature 



|x..(e + i^,r)|cxT-(i-2^"), 



(61) 



as long as nksT ^ 7. The zero temperature scaling 
behavior of the resonance peak (|59p will be masked by a 
finite temperature with a crossover at irksT « 7. These 
scaling and crossover behaviors of the resonance strength 
are features of the coupling of the Majorana edge states 
and the flux qubit4^ 

The finite temperature susceptibility shows a reso- 
nance at e + I/, as shown in Fig. [51 Around the reso- 
nance, the frequency dependence at finite temperature 
will be given by the power law (j60p but with the region 
constrained by TrfesT instead of 7 if tt/cbT > 7. 



VII. HIGHER ORDER CORRELATOR 

So far, we have computed the qubit susceptibilities to 
their first non- vanishing orders and the lowest order self- 
energy correction e 1— > e + iy — ij/2. As a consequence, we 
only used the two-point correlation functions (^(i)^(O)) 
in our evaluations. The next nontrivial corrections to 
the qubit correlators involve the equal position four-point 
correlator of the disorder fields {fj.{ti) iJ.{t2)fJ.{t3) 11(14)) . 



Hence, we confirm that the second order perturbative 
correction is consistent with the the self-energy correction 
calculation. 

The leading correction to the susceptibility Xxx in sec- 
ond order is due to the logarithmic term oc tr^^^logt in 
the correlator (|D33[) and has the form 



J^(2 + v/2)r(|)r(|)e'3"/« 
16£7/4[i(e + 7/2)]i-2A 

^ 1^(7/2)^ + (--^--)' 



(64) 



where we have included the self-energy correction ()54|) . 
and omitted terms without logarithmic divergence. Un- 
fortunately the effects of nontrivial exchange statistics of 
disorder fields are not apparent in this correction. 



VIII. CONCLUSION AND DISCUSSION 

We have proposed a novel scheme to probe the edge 
vortex excitations of chiral Majorana fermion edge states 
realized in superconducting systems utilizing a flux qubit. 
To analyze the coupling we mapped the Hamiltonian of 
the Majorana edge states on the trans verse- field Ising 
model, so that the coupling between the qubit and the 
Majorana edge modes becomes a local operator. In the 
weak coupling regime S s wc have found that the 
ground state expectation values of the qubit spin are 
given by 



gl-2A,,A2A,, 



0, {a' 



3S 

8i< 



(65) 
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Additionally, the susceptibility tensor of the qubit spin 
in the basis x, y, z is given by 

/ 1 -S/e\ 
x{oj) = X..(w) 1 , (66) 
\-6/s {S/e)y 

gj37r/8p(-3-) 

Xxx{uj) = T7-, ; /r,ii-2A A2A ' 

with the real part ly and the imaginary part 7/2 of the 
self-energy given by 

cos^(|)r(|)J^ /o ./^ n 

^ = 2gi-2A,A2A, ' 7/2=(V2-l)j.. (68) 

We see that all of these quantities acquire additional 
anomalous scaling (e/A)^'^'' due to the fact that each 
spin flip of the qubit spin couples to a disorder field /i. 
Similar scaling with temperature appears in interferomet- 
ric setups )^ but using a flux qubit allows to attribute 
its origin to the dynamics of vortices much more easily 
and also gives additional tunability of the strength of the 
coupling. Another effect of the vortex tunneling being 
present is the phase change dcj) = 37r/4 of the susceptibil- 
ity around the resonance4^ This phase shift occurs due 
to the anomalous scaling and the presence of the Abelian 
statistical angle of the disorder field, in view of the fact 
that Xxx is just a correlator of two disorder fields in the 
frequency domain. 

The long wavelength theory which we used is only ap- 
plicable when all of the energy scales are much smaller 
than the cutoff energy of the Majorana modes. This is 
an important constraint for the flux qubit coupled to the 
Majorana edge states. In systems where the time-reversal 
symmetry is broken in the bulk (unlike for topological 
insulator-based proposals^), the velocity of the Majo- 
rana edge states can be estimated to be vm cx vp^/Ep 
and the dispersion stays approximately linear all the way 
up to A. The cutoff of the Majorana modes is related to 
the energy scale of the Ising model A = A H> J. Equat- 
ing J = A and vm ~ 2 Ja, we obtain the lattice constant 
of the Ising model a = vp/Ep = \p, with \p the Fermi 
wavelength. The Fermi wavelength is typically smaller 
than any other length scale, and so the long wavelength 
approximation we have used is well-justifled. For a typ- 
ical flux qubit the tunneling strength 5 is indeed much 
smaller than the superconducting gap, the level splitting 
£ may vary from zero to quantities much larger than the 
superconducting gap. 

Our proposal provides a way to measure properties of 
the non-Abelian edge vortex excitations different from 
the conventional detection scheme that requires fusing 
vortices into fermion excitations. However, none of our 
results for the single fiux qubit can be directly connected 
to the non-Abelian statistics of the quasiparticles, even 
after including higher-order corrections. Thus, it is of 
interest for future research to investigate a system where 
the edge vortex excitations are coupled to two qubits 
such that braiding of vortex excitations can be probedi^ 



Another feature of systems with several qubits worth to 
investigate is the ability of the Majorana edge modes to 
mediate entanglement between different flux qubits. 
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Appendix A: Flux qubit 

The flux qubit which we consider consists of a su- 
perconducting ring interrupted by a Josephson junc- 
tion which is parameterized by its critical current /c, 
its capacitance C, and the self-inductance L of the ring 
threaded by a magnetic flux The Hamiltonian in the 
phase basis reads^i 

H = -^Ec^ + Ej{l - cos(/.) + ^(0 - 27r$/$o)', 

(Al) 

where (j) is the phase difference across the Josephson junc- 
tion and $0 = ^/2e is the superconducting fiux quantum. 
We have introduced the charging energy Ec = e^/2C, 
the Josephson energy E,j = $o-^c/27r, and the inductive 
energy El = 3>o/47r^L. 

The potential energy is given by the last two terms 
of the Hamiltonian (|A1|) . Neglecting for a moment the 
inductive energy, the cosine potential favors states with 
(j) = 27rZ. The transition between these state involves a 
change of the phase difference by 2ti which corresponds 
to driving a vortex in or out of the superconducting loop. 
The inductive energy breaks the degeneracy of the states 
with a different number of vortices n in the loop by fa- 
voring states with n$o ~ ^- When the flux $ is tuned 
close to $0/2, the system becomes frustrated since the 
states (/) = and (j) = 2t: are then nearly degenerate in 
energy. When the inductive energy is smaller than the 
Josephson energy but still large enough such that states 
with more vortices in the superconducting loop are not 
accessible, the potential takes the form of a double-well 
with the minima close to and 27r. These requirements 
are met when El ~ Ej/2Tr^. 

The charging energy Ec describes the influence of 
quantum dynamics. If the level spacing D, ~ ^/SEcEj 
in each well is large enough and additionally the two 
wells are well separated, only the lowest energy states 
|0) and \2tt), which are localized near the classical min- 
ima (/) = 0, 27r, are relevant. Hence, the low energy 
Hamiltonian of the system reduces to Eq. For 
El ^ Ej, the energy detuning of the two minima is 
given by e = 47r^i?i(i — ^/^o) which can be tuned via 
the fiux $ in the superconducting loop. The tunneling 
amplitude is given by (5 cx ex-p{—y^8Ej/Ec)'^ 
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Let us now discuss the experimental parameters for 
the flux qubit. Assuming that the superconductor order 
parameter A is about 1 K, the corresponding coherence 
length is of the order of ^ ^ 1 /im. To avoid the mixing of 
the Majorana edge states, the width of Josephson junc- 
tion needs to be larger than the coherence length which 
is in the range of micrometers. This is consistent with 
most experiments 1^1^ Although the design of the flux 
qubit in Ref. is more complicated than the simplest 
design discussed here, the idea of a 27r phase shift for a 
full vortex tunneling through the Josephson junction is 
the same. Thus, as a concrete example, we quote the ex- 
perimentally achieved parameters from Ref. IH: Ej « 9 
GHz, Ec « 2.5 GHz and El « 0.52 GHz. The tunnel- 
ing amplitude is measured and estimated to be 5 « 369 
MHz 4^ Moreover, the level spacing is estimated to be 
n = 13.4 GHz. 



Appendix B: Effective two-level system 

The Euclidean action of the superconductor phase cor- 
responding to the Hamiltonian (jAip reads 



T/2 



dT 



T/2 



1 1 



2 8Ec 



0' + Vi^) 



(Bl) 



As discussed in App. \^ the double-well potential V{(f)) 
has two energy minima located at (/> = and 27r such that 
V{(j) = 0) = -e/2 - n/2 and V{<j) = 2tt) = +e/2 - n/2. 
For later convenience, we have shifted the potential en- 
ergy by if2. In our discussion, we will assume that the 
level spacing ft = ^/SEjEc is the same at both wells and 
that the potential profile connecting two minima can be 
approximated by V{(f>) ~ Ej{l — cos 0) — ^il. The con- 
crete form of the potential does not affect the qualitative 
feature of our discussion.— 

The action of the Majorana fermions can be inferred 
from the Hamiltonian ([5]) as 

S4, = / —— [iJudiJu + i'ddiPd] , (B2) 

where d = {dr—idx)/2 and B = {dr+idx)/2. The action 
describing the coupling between the phase field and the 
Majorana fermions is given by 



rT/2 ^ 
S^^tj, = i I dr- 

T/2 



dxpe{x,T), (B3) 



where pe{x,T) = i/;'^ {x,t)iP{x,t) with iI;{x,t) — {ipu + 
iipd)/2y/TT is the fcrmion density of the Majorana 
fermions. The origin of this coupling is the electro- 
static energy VQ where V ^ 4>/2 is the voltage from the 
Josephson relation and Q = J^'^ dxpe{x, r) is the charge 
of the superconductor island at one side of the Josephson 
junctioui^ Here we have chosen a gauge such that the 
superconductor phase at the other side of the Josephson 



junction is fixed. Observe that the equation of motion 
of the phase field is not affected by the coupling term 
()B3|) when the integration of the fermion density yields 
no explicit time dependence. 

The total action of the system thus becomes S = + 
S.^, + S'^,.0. From the Euclidean (imaginary time) version 
of Fcynman's path integral, the transition rate reads 



(B4) 



{cj>f\e-"'\^,)=Af [di'] m 



where H is the corresponding Hamiltonian, \4>i.f) repre- 
sent the initial and final phase eigcnstates, and J\f is the 
normalization constant. Because the leading contribu- 
tion to Eq. (jB4p at large times T — >■ 00 comes from the 
lowest-lying energy eigenstates, the Hamiltonian at the 
left hand side can be approximated by an effective Hamil- 
tonian that contains only a few low energy statesi^ 

For the double well potential V{(j)), there exist two 
low energy states that are localized at the two classical 
minima at (f> = and 2tt. By considering the transition 
rates within and between two minima 



(B5) 



for (f>ij ~ 0, 2tt, we would like to show that the effective 
Hamiltonian is a two-level system coupled to the Majo- 
rana fermions. 

To compute transition rates for (f>i = i/iy, we first ob- 
serve that the phase field is mostly localized at one of 
the wells and behaves as a simple harmonic oscillator. 
Therefore, the main contributions to the transition rates 
(|B5[) come from the phase field in the localized states and 
are given by 



^0,0 



and i?o 



-eT/2 



(B6) 



for — (j)f — and 27r states, respectively)^ Notice 
that the phase field and Majorana fermions are effectively 
decoupled when the phase field is localized. 

The other contributions to the transition rates come 
from trajectories of the phase field that contain tunnel- 
ing events between two wells. These tunneling events are 
so called instantons and anti-instantons that occur in a 
very short time interval At l/fi. In the dilute gas ap- 
proximation, each instanton or anti-instanton event cen- 
tered at time Ti contributes to the transition rates with 
a factor 



K±{n] 



■P±{n), P±{t, 



(B7) 



Here, 5 ^ /^c jg i^j-^g action from the tunneling of 

the phase field through the barrier, and P±{t) is due to 
the coupling (jB3p with the approximation that the time 
interval of instanton Ar is small such that the density 
field can be replaced by Pe{x^ Ti). Since the integration of 
fermion density is an integer, we define P+{t) = -P|(t) = 
P_(r) EE P(r). We also note that P^(t) = 1. Moreover, 
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we are not free to distribute the instantons and anti- 
instantons. They have to be alternated in time and the 
first tunnehng event is determined by the initial state.— 
Let us consider the transition rate for Rq.o- Because 
the phase field needs to tunnel an even number of times 
in order to be back to the initial well, the number of 
(anti-)instantons has to be even for a non- vanishing con- 
tribution. A trajectory of the phase field that contains 2n 
(anti-)instantons ordered in time, T/2 > T2n > T2„_i > 
gives the contribution to -Rq.o 



2n 
0,0 



e'^(^JUd^^|e^-'^'"^^^P{n)^ , (B8) 



which is integrated over the centers of (anti-)instantons 
Tj. Here {■■■)ip is the path integral summation over 
fermion fields such that 



{0{i'u,ipd))4, = / [dtp] 0{i!u,il^d) e 



(B9) 



for an arbitrary fermion field combination 0{ipmijJd)- 
The total transition rate can be written as 



-Ro,o — ^ -Rofo' 



(BIO) 



where i?Q g is defined in Eq. (jB6[) . The transition rate 
of i?27r,27r Can be derived in the same manner and takes 
the same form as i?o.o hi Eq. (jBlOp by summing over 

The transition rates between two wells, i?27r,o and 
Ra.2-K, can also be computed by properly counting the 
(anti-)instanton events. The crucial difference is now 
that an odd number of tunnehng events are needed for 
the final state to be in a different well than the initial 
state. The total transition rate of i?27r,o then reads 



j=o 



(Bll) 



by summing over contributions from trajectories with 
odd tunneling events 

Rl::o'-\l^/^j'l[dnleS-'y''^^^P{n)^ . (B12) 

Finally, i?o,27r also takes the same form as R2tt,o with a 
substitution of rI]:^^ = ^2(^0^ [£ ^ -e] in Eq. (|BT21) . 

By using the interaction picture, we can explicitly show 
that the effective Hamiltonian 



(B13) 



reproduces the transition rates within and between two 
wells in Eqs. (jBlOl IBlip up to an overall constant, 
see Eq. Here, Hmf is defined in Eq. ([5]) and 

^x,z Pauli matrix acting on the two- level basis \(f)), 



4> G {0, 2tt}, of the superconducting phase difference 4>'^ 
Note that enters in the Hamiltonian (|B13p together 
with P. This is a consequence of gauge invariance — 
whenever the superconducting phase difference changes 
by 2tt the phase of the Majorana to the left of the junc- 
tion has to be changed by tt. Every physical observable 
has to be gauge invariant that is why (and Ty for that 
matter) always have to occur together with P. In this 
spirit, we define the (observable) qubit degrees of freedom 
as 



(B14) 



see We thus conclude that the two- level Hamilto- 

nian (jBlSp together with the identification (jB14p gives 
the effective low energy description of the system. 



Appendix C: Correlation functions of disorder fields 

The one-dimcnsional critical transverse-field Ising 
model is a conformal field theory (CFT) with central 
charge c — 1/2. This CFT contains the following pri- 
mary fields: 1, e = itpip, s, and ^. Here 1 is the identity 
operator, e is the energy field (a product of the right and 
left moving Majorana fermion fields ip and "0), and s is 
the Ising spin field with its dual field The dual 

field fi is also called the disorder field and has the same 
scaling behavior as the Ising spin field s at the critical 
point. On the lattice, the disorder fields fj, are non-linear 
combinations of Ising spin fields s and reside on the bonds 
of lattice Ising model. They are hence not independent 
of the Ising spin field s. 

In the continuum and in imaginary time, the two-point 
correlation function of disorder fields can be obtained 
from CFTil 



1 



[{Zl - Z2){Z1 - Z2)] 



(CI) 



with Zi = Ti + ixi and Zi = Ti — ixi . 

Following Ref. [13, the real-time correlators can be ob- 
tained by analytical continuation r £,+it. Here ^ 0+ 
is introduced to ensure the correct phase counting and is 
important for the Abelian part of the statistics. The 
equal position two-point correlation function is given by 



{fl{ti,Xo)fl{t2,Xo)) = 



By using the identity 



lim 



1 



i^ + i{ti-t2)) 



2A„ 



-i sgn(*)7r/8 



5^0+ (C + zt)l/4 |t|l/4 



(C2) 



(C3) 



one obtains the two-point correlation function in the form 
of Eq. (El]). 

The four-point correlation function of /i's can be ob- 
tained in a similar manner. In imaginary time, the cor- 
relation function is given hy^ 
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{H{zi,zi)fl{z2, Z2)H{Z3, Z3)^(Z4, Zi))" 







-1x1 + 


212^34^14223 


V 2 



(C4) 



where x = (212234/2:13224) is the conformally invariant cross ratio, and the absolute values should be understood as 
l-^yl" ~ i^ij^ij)"^'^ ■ Because we are interested in tunneling at a single point, we can set Xi = 0. In this limit the 
four-point correlation function can be evaluated to be 



(Ai(zi)^(z2)/i(23)/i(z4))^ 





213224 


1/2 












212234214223 












< 


213224 


1/2 

|1 


- X 




214223 




212234214223 


= 


212234213224 




213224 


1/2 

Ixl 




212234 


1/2 




212234214223 




214223213224 





1/2 



for < X < 1 

for X < (C5) 

for X > 1 



The real-time correlation function can be obtained by first taking a square root of Eq. (|C5[) followed by the analytical 
continuation, + itij^ 



^Fi2iti,t2,t3, t^) [6'(1324) -t- 6'(1423) + 61(2413) 61(2314) + 61(3241) + 61(3142) + 61(4132) + 61(4231)] 
+^13 {h MM. h) [6'(1234) + 61(1432) + 61(2143) + 61(2341) + 61(3214) + 61(3412) -t- 61(4123) + 61(4321)] 
+i^i4(ii, i2, ^3, U) [6'(1243) + 61(1342) -I- 61(2134) + 61(2431) 4- 61(3124) 4- 61(3421) + 6/(4213) + 6/(4312)] , 

where 6{abcd) — 1 for ta > % > tc > td and is otherwise zero. The corresponding functions Fij are given by 

Fl2{tl,t2,t3, ti) = 
Fl3{ti,t2, ^37 ^4) = 

F,iih,t2, h, U) = ^ ^^^^ _ ^^^^^/^^^ ^ ^^^^ _ ^^^^^/^^^ ^ ^^^^ _ ^^^^^/^^^ ^ ^^^^ _ ^^^^^/^ 



(C6) 







[en 


-i(ii 


-t2)V/%^ 


~i{t3 








-z(ti 








-l{t2 


-t3)m^^ 


-l{t2-U)]^/^ 








-i[ti 




-i{t2 








-z(ti 


-i2)]i/4[eH 


- i(ii 




~l{t2 


-i3)]i/4KH 


-z(t3-t4)]l/4 








-i(ti 




hi(t2 


-^3)]^/^ 





(C7) 



Here i^i2, F13, and F14 are the three characteristic 
functions appearing in the fourth-order correlation func- 
tions. For an Abelian state, they usually appear in quasi- 
symmetric combinations and exchanging two of the times 
alters various phase factors, which is a characteristic of 
fractional statistics. For the current non-Abelian case, 
however, exchanging two of the times not only alters 
phase factors but can also change the form of the correla- 
tion function from one of the characteristic functions to 
another. This is a special feature of non-Abelian statis- 



Appendix D: Second order correction to {a^{t)a^{0))c 

Because our ultimate goal is to compute the qubit sus- 
ceptibility, we are interested in the correlator with t > 
in the long-time limit < — > cxj. Let us first recall the 
perturbative part of Hamiltonian (j27p in the interaction 



picture: 

V{t,xo) = -^-mr^it)+r-{t)]. (Dl) 

Since the vortex tunneling in or out of the superconduct- 
ing ring directly couples to the disorder field of the Ising 
model cr^(i) = ii{t)T^{t) in the transformed basis, the 
evaluation of the second order correction for the correla- 
tor (cr^(t)(T^(0)) requires the knowledge of the four-point 
correlation function derived in Appendix [C] 

We expand the 5 and 5^-matrices in ([5T|) to second 
order with insertions at times ti and ^2- Nonzero contri- 
butions to the correlator come from three regions: (A) 
t > > ti > t2, {B) t > h > > t2 and (C) 
t > ti > t2 > 0. These three regions are shown in Fig. |31 
In what follows, we will evaluate the second order con- 
tributions from each region in the long-time limit. 
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FIG. 3: The integral domains for regions A, B and C in the 
ti and t2 coordinates used in Appendix ID] 



1. Region A: t > > U > t2 

The contribution from region A is given by 

{a^t)a^O))f = {-if f dh r dhiA, (D2) 

J — oo J —oo 

with the integrand 

Ia^ + {a^t)a^0)ViV2)o + (^2 Vi (7^(0^7^ (0))o 



(D3) 



where Vi = V{ti) is a shorthand notation. The plus and 
minus signs come from the location of the insertions. The 
plus sign corresponds to having both insertions located 
on the same branch (either forward S or backward S^) 
while the minus sign corresponds to the situation where 
the two insertions are located on different branches. 

Because only certain ordcrings of insertions of raising 
and lowering operators or t~ , coming both from the 
interaction term (jPlj) and the , give non- vanishing con- 
tributions, the integrand is given by 



1 1 /. 



-iet ie{t2—ti) 



{^i{t)^i(S))^i{h)^i{t2)) 



+ e-^^*e'^(*i-*^)(M(t2)M(ii)M(i)M(0)) (D4) 



^+ietgie(ti-t2) 
,+iet ie(t2-ti) 



{ll{t2)ll{t)ll{Q)ll{h)) 

{^i{h)^i{t)^i{Q)^i{t2)). 



Here, the four-point correlation function can be read off 
from Eq. (jC6|) and simplified using the identity (|C3I) . 
Remarkably, these correlators have the same time de- 
pendence function and differ only by phase factors. This 
feature is characteristic also to regions B and C. After 
some algebra, the integrand simplifies to 



To evaluate the integral (jD2p . we first simplify it by introducing new variables such that <i = —tT and t2 = —t{T+T) 
with the new integrating domain < r < cxj and < T < 00. The second order correction from region A becomes 

where we have introduced a regularization factor exp(e77ti), with 77 — > 0^. 

The integral in Eq. (jD6[) will not generate any oscillatory dependence but is divergent when both T and r are large. 
It is thus convenient to separate the algebraic part of the integrand into three parts 

(l + T)V4(T + r)i/4 1 ^3/4 1 ^3/4 



(l + r + T)i/4(TT)i/4 ^1/4 4(T-f r)(l + r)' ri/4' ^-^ 4{T + t){1 + T)' 

Combined with the exponential prefactor, the integration of Iai is regular, the integral of IA2 diverges linearly while 
that of Ias diverges logarithmically. 

Integrating /^^ with all the exponential prefactors gives 

-drdT = ^ , 'J ex --L-[-zik2 + + ^ 0+ (D8) 



ri/4 2?7(et)7/4(j + ^)3/4 (et)7/4 ^ 2r/ 

Since the the linear long time divergence is purely imaginary, it does not contribute to the correlation function. 
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In the long-time limit, the integrals of Iai and with all the exponential prefactors can be carried out to the 
lowest order in l/{st) and are given by 



r(|)(7r-2(l + log(8))) 



r(|) (31og(8rf) - (3^^2)6-^ +37-4) 



12(£t)V4 



(D9) 



(DIO) 



We now add the real parts of the three integrals (|D8|) , (|D9|) , and (jPlOp and then multiply them with the prefactors 
in (jD6p . The result is the leading long-time contribution from region A to the qubit spin correlator: 

i^^^'itWmf |r(|)(7 + 3.-3,-181og(2)-31og(st))| _ ^^^^^ 
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In the long-time limit, the leading contribution is given by t he term oc t-V4 iog(,t). 



2. Region B: t > ti > > t2 

The contribution from the region B is given by 

^(Oa-W)*,'^ = i-if f dh f dt2lB, (D12) 
Ja J -00 

with the integrand 
Ib= + {<j''{t)V,<j-{0)V2)o + {V2Via-{t)a^0))„ 



(D13) 



After ordering the the raising and lowering operators 
or and using Eq. (jC6p . the integrand reads 



Ib =e~'^l^lB, + e+'^/^/s, - 11^ - , (D14a) 



where the two integrand functions are given by 



' "4(t - ti)l/4(i _ t2)l/4(t^)l/4(_i2)l/4 ' 
^ "4(i-ii)l/4(t_t2)l/4(t^)l/4(„i2)l/4' 



(D14b) 
(D14c) 



with X* denoting complex conjugate of x. Again, the 
four-point correlators of ju's in region B have the same 
functional form up to phase factors. 

To evaluate the integral of Ibi , we introduce new vari- 
ables xi and X2 with ti ~ t{l — xi) and t2 ~ -~tx2 such 
that 



We can then split the integral Bi into an oscillatory contribution i?f and a non-oscillatory one B^'-^ . 

Since the non-oscillatory contribution from (|D15|) is dominated by xi ^ X2 ~ 0, we can expand the integrand 
around this point to get the leading contribution. Because we are interested in the correlator in the long-time limit, 
we then deform the integration contour in the complex plane such that both xi and X2 change from to — ioo. The 
leading non-oscillatory contribution is given by 

MO ^ / dx, / dx2e-^^'*(^^+^^' ( — Attt + ^ ' ^' 1 = — Ht^ ( -1 + — W ) • (D16) 

' 4 Jo Jo \{xiX2Y'^ 4 ; 4£3/2 1^ MeH-^J ^ ^ 

The oscillatory contribution Bi is dominated by xi w 1 and X2 ~ 0, we can thus expand the integrand around this 
point to get the leading contribution. Again we are interested in the correlator in the long-time limit and thus deform 
the integration contour such that xi varies from 1 — ioo to 1 and X2 varies from to —ioo. After these transformations 
i^P evaluates to 
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where ui = xi — 1. 

Summing up, the leading contributions to Bi are 



4e3/2 



64e2t2 



+ e 



-.etCos(f)r(|)r(|)r(|) 



v^r(f)(rf)i/4 



(D18) 



The leading non-oscillatory contribution of Bi is a constant while the leading oscillatory contribution has a power 
law decay oc t~^^^. 

To integrate /s^, we again use the variables ti ~ t{l — xi) and t2 = —tx2 such that 



(l-Xi+X2)l/4 



(:ri)l/4(l + :r2)i/4(i_ 2.^)1/4(2.^)1/4 



(D19) 



Once again, the non-oscillatory contribution is dominated by xi ^ X2 ~ 0. We expand the algebraic part of the 
integrand around xi — X2 = 0, deform the integration contour such that xi runs from to —ioo and X2 from to 
ioo, and get 



B. 



NO 



3/4 3/4 ■ 



axi I ax2e <^ -h ^ j> ^^^^^ 



-/4 1/4 

I •^'2 



(D20) 



To evaluate the oscillatory part B2 of B2, we expand the integrand around Xi — \ and 2:2 = for the leading 
contribution. The necessary deformation of the integration contour is now given by Xi changing from 1 — ioo to 1 and 
X2 from to ioo. The leading oscillatory contribution from Eq. (|D19[) is now given by 



^2^3/2 _ /.O 

" ' dui / dx2e 



with ui = a;i — 1. 

The final expression for B2 is 



J (x2-^l)^/^(-^^l)^/^ 
I 1/4 



re * e , 



4eV4 y^ii/4 



(D21) 



B2 



4e3/2 



64e2t2 ^ ^(et)i/4 



(D22) 



Similarly to _Bi, the leading non-oscillatory contribution of B2 is a constant, while the leading oscillatory contribution 
has a power law decay '-^ t~^^^. 

From Eq. ()D12I) and (jD14l) , the leading contributions to the qubit spin correlation function from region B is given 

by 



2e3/ 



4/ M 9*sin(V4) , 3cos( f)r(|) ^_.(,,+.;4)) , r(|)eW%^_,,, 



1 — cos(7r/4) 



64£2i2 8y2r(^)(rf)i/4 



^/^(e<)'/ 



— (e-*-e-') 



(D23) 



3. Region C: t > h > tn > 

The integral in region C reads 

(a^(t)a^(0))g^ = i-tf fdh r dt2lc- 
Jo Jo 



(D24) 



We calculate the integrand Ic in a similar way to regions 
A and B. We get 



(D25a) 



with the two integrand functions being 

J2g-»etg»e(ti-t2)(i _ ^2)1/4(^^)1/4 



4(t-ti)V4(t)l/4(i2)l/4(i^_i2)l/4' 
j2g»etg-»s(ti+t2)(^_^^)l/4(^^)l/4 
^4(<-ti)l/4(t)l/4(i2)l/4(i^_f2)l/4- 



(D25b) 
(D25c) 



To integrate Ici , we make the variable transformation: 
h = <(r+l/2 + T/2) andt2 = <(r+l/2-T/2). In terms 
of the new variables, the integral of Ic\ reads 
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_5^f^^e^ f\ e-*- ^i/2-./2 (^^2 _ y ^ r/2) V4(l/2 + T + r/2) 



The integration over T can be carried out exactly with the result 

= i F(irX '"wi^i-^'^^^(^-i'2'4'(,TT7j J' 

where 2-^1(0;, /3; 7; x) is the Gaussian hypergeometric functioni^ 

We deform the integration contour in Eq. (|D27p such that r goes from to +ioo and then back from 1 + ioo to 1. 
The leading contribution in the long-time limit is dominated by the region near the real axis. The expansion around 
X = leads to an oscillatory contribution while the expansion around x = 1 leads to a non-oscillatory contribution. 
To the lowest few orders, the asymptotic behavior in the long-time limit is given by 

^1^^ 



4£3/2 \^ 64£2t2 



e3W8r(|)i3/4 e^*''/«r(|)(61og(rt)- (6 + 3i)7r + 67- 14 + 361og(2))l ^^^^-^ 



r3/4 



12e7/4ii/4 J 



The oscillatory contribution contains a power-law divergent t^^^ term. As we discuss later, this term contributes to 
the shift of the resonant frequency and to the damping for the (cr+(t)cr~(0)) correlation function. 

To integrate /cs, we first change the integration variables to r and T defined by ti — t(T + r/2) and t2 — t{T — t/2) 
such that the integral separates into two parts 

^2^3/2g.st j- ^1/2^^ ^2T^^ ^1 /■2-2T^_] f^-aetT {I - T + T /2f'HT + T /2)^/^ 



After changing T — > 1 — T in the second integral and then introducing X = 2T , this equation simplifies to 

^ -^""^ Y Jo Jo '"(2-X-r)V4(X-r)V4(,)i/4j ' (^^O) 

Again, we deform the integration contour in the integral over X with X changing from to — ioo and then from 
1 — ioo to 1. Now the oscillatory contribution comes from X ^ Q while the non-oscillatory one from X ^ 1. By 
expanding the integrand around these two points, we get the leading contributions: 

9 + e^^-^e— /«) (-1, f ; |; -l) r(|) \ 

The Gaussian hypergeometric function evaluates to 2F1 {—j, |; f ; — l) ~ 1.102. 

From Eqs. (|D24p and (|D25p . we obtain contribution to the qubit correlation function from the region C: 

(a-(i)a-(0))[?) = - (e— /4 + _ c^) 

S'm' f. , 9^sin(f)^ , r_;se-*\ 5^V{1) 



j2g-«£tg-^7r/8(-g~»7r/4 _^ w r(I)(6 log(rf) - (6 + 3i)7r + 67- 14 + 36 log(2)) ^ (D32) 
4^W74 \ 12 J 

^'r(|)2(e-W4 + 1) (e»£«e»V8 ^ ^-^et^-^./B^ ^p^ (_ l ^ 3. |. p(7^ 

4eV4ti/4 



4. Final result for (cr^(t)cr^(0))i^-' butions from all the three regions, given by Eqs. (|Dlip . 



The second order correction to the correlation function 

f 21 

(fT^(<)cr^(0))c can be obtained by adding up the contri- 
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(jD23[) and (|D32[) and then subtracting (cr^)^ as calcu- in the long-time hmit is 
lated in Eq. PO)) . The full expression for the correlator 



_ ^^r(|)e-/«e---^ r ^-,W4. (7 + 3.-37-181og(2)) / . 



'^'r(|)2r(|)e-'^* /2e^ 3cos(f)r(|)e-W4 2 i^i (-j |; |; -l) (e"^'"/' + e-""/«) 



4e7/4ti/4 \^ ^ 472r(f )r(|) 

'^^r(|)V^* / [7 + 37r-37-181og(2)]e— /8 3cos(f)r(|) 
4£W/4 I 8r(f) ^ 4\/2r(J^) 

2r(|)e-/8 2cos(f)2Fi(-i,|;|;-l)r(|) 



(D33) 



This result agrees well with numerical evaluation of the 
integral. A power law divergence ^ t^^'^ and a logarith- 
mic contribution ~ \og{et)/t''-^'^ dominate the long-time 
behavior of the correlator. However, this logarithmic con- 
tribution will be cut off either by the induced damping 
or by a finite temperature. 

A heuristic way to see that the term diverging as t^^^ 
corresponds to self-energy correction is to add it to the 
zeroth order correlator of (cr^(i)o'^(O)) given by (|4T|) . The 
sum of these two terms equals to 



-iTz/8'_ 



1 



^(2cos2(f)-i^)t 



(D34) 



, 1 1 ~ — i — )t , 



with ly and 7 the same as in Eq. (|5ip . It then becomes 
apparent that (jD34[) is exactly the expansion of the renor- 
malized correlator (|55p to the second order in S 



in/S 



tl/4 



-i(e+y)t-7t/2 



(D35) 



We thus conclude that the explicit evaluation of the 
higher order correction gives a result consistent with the 



self-energy calculation. 



Comments on leading contributions of higher 
orders 



The leading contribution to the second order correc- 
tions comes from region C with integration of Ci when 
the integration variable r is around r = 0, cf. Eq. ()D26I) . 
Since t = {ti — 12) 1 1, this expansion to the zeroth order 
is equivalent to making an operator product expansion of 
fj,(ti) ii{t2) for ti « t2 in the four-point correlation func- 
tion of /i operators4i In the n*'' order of perturbation 
theory with insertion times ti, . . . , we expect that the 
most divergent contribution arises when all the insertion 
times belong to the interval [0,t]. By ordering the times 
ti > t2 > ■ ■ ■ > tn and using the operator product expan- 
sion for the pairs A't2i-iA*t2j for i = 1, . . . , n/2, we get a 
perturbative structure resembling Wick's theorem. The 
resummation of these terms would give the contributions 
for the self energy which we calculated in Sec. IVI Al 
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The basis for the two-level system, \(f) = 0) and \cj> — 2tv), 
used in the App. [B] directly denotes the superconductor 
phase difference across the junction. This basis is in turn 
different from the basis |0) and |27r) defined in the main 
text, right above Eq. 0. They are however related by 
|0) = = 0) and \2tv) = P\(l> = 2tv). 



